## This code simply: 
## (1) Typesets some paragraphs with descriptive data
## (2) Produces the graph in Appendix I that describes the data
## Requires: dataset-Nbydatemun.RData -> a summary by date of cummulative mcmv units and municipalities
## Saves figure in a "Figures" subfolder, which is assumed to exist

rm(list=ls())
library(plyr)

load("Data/dataset-Nbydatemun.RData")
tmp <- mcmv.summary

## (1) Typeset paragraphs that describe the delivery of units over time 
## This paragraph is included in the DiD section of the main body of the paper

bfirst <- tmp$date[1]
b2010 <- tmp[max(which(tmp$date<="2010-12-31")),]
b2014 <- tmp[max(which(tmp$date<="2014-12-31")),]

t2015 <- tmp$cumunits[max(which(tmp$date<="2016-01-01"))]-tmp$cumunits[max(which(tmp$date<="2015-01-01"))]
t2016 <- tmp$cumunits[max(which(tmp$date<="2017-01-01"))]-tmp$cumunits[max(which(tmp$date<="2016-01-01"))]

t2020 <- tmp$cumunits[max(which(tmp$date<="2021-01-01"))]-tmp$cumunits[max(which(tmp$date<="2020-01-01"))]
t2019 <- tmp$cumunits[max(which(tmp$date<="2020-01-01"))]-tmp$cumunits[max(which(tmp$date<="2019-01-01"))]

b2020 <- tmp[max(which(tmp$date<"2021-01-01")),]

cat("The first beneficiaries received their units in late 2009, and by the end of 2010, ",
    round(b2010$cumunits/1000)*1000,
    " had been delivered in about ",
    round(b2010$cummun/100)*100,
    " municipalities. The program expanded in the next cycle, reaching more than ",
    round(b2014$cumunits/1000)*1000,
    " beneficiaries in ",
    round(b2014$cummun/100)*100,
    " municipalities. The MCMV delivered ",
    floor(min(t2015,t2016)/1000)*1000, 
    " units more in both 2015 and 2016, but lost momentum after that, with the acute political and economic crisis that engulfed the country. By the end of 2020 ",
    round(b2020[,"cumunits"]/1000000,1),
    " million units had been delivered in ",
    b2020[,"cummun"],
    " municipalities in the country.")


# (2) Figure in Appendix I describing delivery of units over time
pdf(file="Figures/fig-I1.pdf",
    width=8,height=6)
par(mar=c(4,4,2,4))
plot(as.Date(tmp$date),tmp$cumunits,type="l",xlab="",
     ylab="Cummulative Units Delivered  (Faixa 1)",
     bty="n",xlim=c(as.Date("2009-01-01"),max(tmp$date)))
text(as.Date("2009-09-01"),1400000,labels="Pres. Cycle 1",pos=3,xpd=NA)
text(as.Date("2012-10-01"),1400000,labels="Pres. Cycle 2",pos=3,xpd=NA)
text(as.Date("2016-10-01"),1400000,labels="Pres. Cycle 3",pos=3,xpd=NA)
par(new=T)
plot(as.Date(tmp$date),tmp$cummun,type="l",
     xaxt="n",yaxt="n",bty="n",xaxt="n",xlab="",ylab="",
     lty=2,xlim=c(as.Date("2009-01-01"),max(tmp$date)))
axis(side = 4, at = pretty(range(tmp$cummun)))      # Add second axis
mtext("Cummulative Municipalities Treated (Faixa 1)", side = 4, line = 2)       
abline(v=as.Date("2010-10-03"),col=gray(0.6),lwd=1)
abline(v=as.Date("2014-10-05"),col=gray(0.6),lwd=1)
abline(v=as.Date("2018-10-07"),col=gray(0.6),lwd=1)
legend(x="bottom",horiz=T,legend=c("Units","Municipalities")
       ,bty="n",lty=c(1,2),inset=-0.15,xpd=NA)
dev.off()

